Thermostability modification of β-mannanase from Aspergillus niger via flexibility modification engineering

Introduction β-Mannanases can hydrolyze mannans, which are widely available in nature. However, the optimum temperature of most β-mannanases is too low to be directly utilized in industry. Methods To further improve the thermostability of Anman (mannanase from Aspergillus niger CBS513.88), B-factor and Gibbs unfolding free energy change were used to modify the flexible of Anman, and then combined with multiple sequence alignment and consensus mutation to generate an excellent mutant. At last, we analyzed the intermolecular forces between Anman and the mutant by molecular dynamics simulation. Results The thermostability of combined mutant mut5 (E15C/S65P/A84P/A195P/T298P) was increased by 70% than the wild-type Amman at 70°C, and the melting temperature (Tm) and half-life (t1/2) values were increased by 2°C and 7.8-folds, respectively. Molecular dynamics simulation showed reduced flexibility and additional chemical bonds in the region near the mutation site. Discussion These results indicate that we obtained a Anman mutant that is more suitable for industrial application, and they also confirm that a combination of rational and semi-rational techniques is helpful for screening mutant sites.


Introduction
Plant cell walls consist of hemicellulose and cellulose, which are both abundant polysaccharides in nature (Katrolia et al., 2012). Mannans are the main components of hemicelluloses, are widely found in plant seeds and endosperm, and are formed by mannose linked via ß-1,4 glycosidic bonds (Singh et al., 2018;von Freiesleben et al., 2019). Mannans can be classified into pure mannans, galactomannans, glucomannans, and galactoglucomannans based on their composition (Wang et al., 2015). Locust bean gum (LBG) and guar gum are commonly mixed substrates for mannanase, and the ratio of pure mannose to galactose in LBG and guar gum is 4:1 and 2:1, respectively (Cao et al., 2018).
Endo-ß-mannanase (EC 3.2.1.78) is present in plants, animals, and microorganisms (Yamaura and Matsumoto, 1993;Zhou et al., 2014;Leerawatthanakun et al., 2022). It can be classified into glycoside hydrolase (GH) families 5, 26, and 113 based on its sequence (Lombard et al., 2014). 1 But in recent years, it is also considered to be a part of the GH134 family Sun et al., 2021). Mannanase degrades mannan to mannooligosaccharide (MOS) by cleaving the mannan backbone, and the MOS promotes the growth of probiotic bacteria such as Bifidobacterium adolescentis in the mammalian gut (Lv et al., 2013;Kumar Suryawanshi and Kango, 2021). In addition, mannanase is also applied in the fields of dye decolorization, pulp bleaching, detergent enhancement, coffee viscosity reduction, and oil drilling (Chauhan et al., 2012;Li et al., 2014;Sakai et al., 2017;Wang et al., 2018;Guo et al., 2021). However, most of the above industrial applications are currently carried out in high-temperature environments, thus, heatresistant mannanase is the preferred enzyme in the feed granulation process, saccharification, and fermentation processes (Turner et al., 2007;Bhalla et al., 2013). However, the optimum temperature for microbial-derived mannanase is generally 40-60r, thus, obtaining a strain of ß-mannanase with excellent thermostability can create significant benefits.
Directed evolution is a traditional technique for accelerating the evolution of enzyme molecules by mimicking natural environmental changes, but random mutations in genes tend to produce large libraries of mutants; therefore, obtaining an excellent mutation by directed evolution is time-consuming and costly (Qiu et al., 2021). Meanwhile, many rational methods for improving the thermostability of enzymes have been reported, such as the design of surface amino acid charges, N-or C-terminal enzyme modification, mutation-free energy calculation, and molecular dynamics simulation (Mahanta et al., 2015;Srivastava et al., 2016;Herrera-Zúñiga et al., 2019;Liu et al., 2021). Rigid flexibility sites (RFS) are one of the protein rational design strategies. They could efficiently screen for mutation sites and improve enzyme thermostability based on the theory that enzyme inactivation at high temperatures was caused by the unfolding of its flexible regions due to heat (Yu et al., 2017). Thus, the flexible regions must be identified and rigidified to maintain enzyme activity at high temperatures. The B-factor and Gibbs unfolding free energy change ( G) are commonly used screening parameters in the rational design. They reflect the mobility of an atom during nuclear magnetic resonance and the conformational stability of the mutant compared to its wild-type (WT) counterpart, respectively (Kellogg et al., 2011;Yuan et al., 2018). Molecular dynamics simulations, called "computational molecular microscopy, " are a well-established technique that can be applied to study the structural, energetic, and thermodynamic properties of various molecules at any temperature and pressure (Feng et al., 2019;Mazurek et al., 2021).
To the best of our knowledge, this is the first study to achieve heterologous expression of an acidic GH5 family ß-mannanase from Aspergillus niger (accession number: XP_001390707) in Pichia pastoris. Then, B-factor and G were used to screen for five positive single-point mutants, all of which showed a slight improvement in thermostability. Furthermore, the obtained multipoint combinatorial mutant showed better thermostability. Finally, the mutational effects were studied using multiple sequence alignment and molecular dynamics simulation. As a result, the present study described a new acidic ß-mannanase with good heat resistance for industrial applications and showed a successful combination use of rational and semi-rational design in the screening of thermophilic enzymes.

vectors, and materials
The ß-mannanase gene was synthesized by Shanghai Sangon after codon optimization. pPICZ vector (preserved in our laboratory) was used to form a pPICZ-Anman recombinant vector with Anman. Escherichia coli JM109 and P. pastoris X33 were kept in our laboratory for cloning and heterologous expression of mannanase, respectively. Restriction endonuclease and DNA polymerase were purchased from Takara (Dalian, China). PCR purification kit was purchased from GENEWIZ (Suzhou, China). LBG was purchased from Sigma (St. Louis, MO, USA) as substrate. The rest of the chemical reagents not specified were analytically pure.

Homologous modeling and calculation of G
Because of the high protein sequence similarity between 3WH9 (PDB ID) and Anman , we selected 3WH9 as the template for Anman using the SWISS-MODEL server for homology modeling (Supplementary Figure 1). The G of Anman calculated using the Discovery Studio 2019 software (Fu et al., 2017). Briefly, using the DeepView software to complement the missing side chains and residues, the Charmm36 force field was imparted to the protein, and G was calculated under pH-dependent conditions.
G was calculated as shown in Equations (1) and (2): G folding indicates the Gibbs unfolding free energy change of WT or mutant in the folding and unfolding states, and the mutant is predicted to be stable if it possesses negative G mut , and the opposite is predicted to be unstable.

Selection of mutation sites
The B-factor of wild-type mannanase Anman was calculated by ba2r, 2 and we subsequently obtained the secondary structure information of Anman using the PDBsum database. 3 To reduce the loss of enzyme activity caused by the mutation effect, the average distance between catalytic residues and ß-turns in loops was measured using Pymol. Subsequently, we further narrowed down the range of mutation sites according to the amino acid preference of ß-turns and then combined it with G screening to obtain hotspots. Meanwhile, Anman was compared with ß-mannanase with 13 excellent thermostability by multiple sequence alignment to explore semi-rational design's role in finding hotspots.

Enzyme expression and purification
The recombinant vectors of Anman and mutants were linearized by restriction endonuclease Sal al restriction endonuclease. The linear vector was transferred into P. pastoris X33 and incubated in yeast extract peptone dextrose (YPD) solid medium containing 100 mg/mL bleomycin (zeocin) at 30 • C for 72 h. Colonies were picked and incubated in 30 Ml YPD liquid medium at 30 • C and 220 rpm for 24 h, and cells were then collected by centrifugation at 8,000 rpm and transferred to 30 mL YP medium containing 0.5% methanol at 28 • C and 220 rpm for 72 h. Methanol was added every 12 h to maintain its concentration at 0.5% for induction of enzyme production.
The induced fermentation broth was centrifuged at 4 • C at 8,000 rpm, and the collected supernatant was dialyzed for 16 h in 0.1 M NaH 2 PO 4 -Na 2 HPO 4 containing 20 mM NaCl (pH 5.5, buffer A) to remove salt and metal ions. Anion exchange chromatography was used for the purification of crude enzymes. Briefly, after preequilibrating the anion exchange column (Hi Trap TM 1 mL Q HP) with Buffer A, the dialyzed crude enzyme was loaded onto the column, followed by a linear gradient elution using Buffer A with Buffer B (Buffer A with 1 M NaCl), and the collection containing the enzyme activity was utilized for the subsequent study. The purification of the samples was analyzed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) (Laemmli, 1970). The Bradford method was used to determine the protein concentration (Bradford et al., 1976).

Determination of mannanase activity
The determination of ß-mannanase activity was performed by the 3,5-dinitrosalicylic acid (DNS) method. Briefly, the purified enzyme solution and LBG were diluted and dissolved with 0.1 M acetic acidsodium acetate buffer (pH 3.0, Buffer C), and the reaction system consisted of 80 µL of pure enzyme and 160 µL of LBG (0.3 mg/mL), incubated at 37 • C for 30 min and then terminated by adding 200 µL of DNS. After boiling for 5 min, 560 µL of water was added to the reaction system, and the absorbance of the reaction solution at 540 nm was measured using a microplate reader. Under the same conditions, the absorbance of mannan oligosaccharides was used as the standard curve. One unit of ß-mannanase activity was defined as the amount of enzyme that releases 1 µmol of reducing sugar per minute under the above conditions. All experiments were performed in three parallel experiments.

Determination of biochemical properties of mannanase
The specific activity of mannanase at different temperatures (50-80 • C) was determined using the method mentioned above, and the temperature with the highest specific activity was determined as its optimum temperature. Thermostability was determined by incubating mannanase at 70 • C for different times (10-60 min) and then reacting with LBG at 37 • C to determine the residual activity. Tris-HCl (pH 2.0) and acetic acid-sodium acetate (pH 3.0-7.0) were used to determine the optimum pH of mannanase, the specific activity of mannanase was measured at different pH conditions (2.0-7.0), and the pH with the highest specific activity was determined as the optimum pH.
To determine the t 1/2 of mannanase, the enzyme was incubated at 70 • C for different periods (5-90 min) and then reacted with LBG to determine its remaining specific activity. The kinetic parameters were determined using 1-10 mg/mL LBG as the substrate (concentration gradient of 1 mg/mL) with the pure enzyme at 37 • C and pH 3.0 to determine specific activity. The K m and K cat parameters were obtained by non-linear fitting using the Origin 2017 software according to the specific activity and LBG concentration. All the above experiments were performed in three parallel experiments.

Molecular dynamics simulation
Molecular dynamics simulations of wild-type Anman and its mutants were performed using Gromacs 2018. Briefly, the missing atoms were added to the PDB file of Anman using the DeepView software before adding the Amberff14sb force field, followed by adding the TIP3P water model to the ortho-dodecagonal box and neutralizing the charge in the system with a sufficient amount of sodium ions. Next, the entire system was minimized using the fastest descent method (1,000 steps), followed by NVT and NPT simulations for 200 ns, respectively, using the LINCS algorithm to limit the heavy atoms and all bonds. The above processes were performed under periodic boundary conditions (PBC). Finally, the system was simulated at 370 K for 50 ns, with a trajectory saving interval of 2 fs. The parameters of root mean square deviation (RMSD), root mean square fluctuation (RMSF), and hydrogen bonds during the simulation were analyzed by the saved trajectory file, and VMD carried out the analysis of salt bridges.

Determination of melting temperature
The protein melting temperature (T m ) represents the thermodynamic stability and reflects the conformational changes of the enzyme at high temperatures. In this study, T m was determined using the Nano-Differential SCNAning Calorimeter (Nano-DSC, Waters, USA) instrument. Briefly, the pure enzyme was diluted to a concentration of 0.2 mg/mL with Buffer C and loaded into a capillary cell. The scanning temperature range was 20-100 • C with a ramp rate of 1 • C/min, and Buffer C was used as the control. The data were analyzed using the Nano Analyze software.

Constraint network analysis
Constraint network analysis (CNA) simulates the enzymes as body-and-bar networks, where the forces between atoms are replaced by bars. The number of bars represents the strength, and E hb measures the magnitude of the energy contained in the force. CNA simulates a thermal unfolding process with increasing temperature, and the E cut−hb decreases as the temperature increases. During the CNA simulation, there is a sudden decrease in the rigidity order parameter (P ∞ ), which means that the network changes from rigid to flexible, just like the protein folding-unfolding transition, and the E cut−hb is the phase transition point. The temperature of conformational transition (T) can be calculated by using Equation (3): Frontiers in Microbiology 03 frontiersin.org Constraint network analysis was accomplished on a web service. 4 In this study, the range of E hb was set from -0.1 to -0.6 kcal/mol, and the decreasing gradient of E cut−hb was -0.1 kcal/mol.

Determination of flexible regions
To identify the more flexible regions in Anman, the average B-factor of α-helices, β-sheets, and loops was calculated, which were 20.5, 17.7, and 21.4, respectively (Supplementary Figure 2). The regions with a higher B-factor level tended to contain more flexible residues; therefore, the present study focused on the mutation range in loop regions, and it has been also confirmed that the modification in loop regions could improve thermostability (Cheng et al., 2012). Due to the different flexibilities of the enzyme surface and internal loops, the B-factors of residues exposed on the surface vs. those buried in the interior were calculated as 24.1 and 18.2, respectively. Meanwhile, the calculated average B-factor of residues within 5 Å of the catalytic residue was only 14.8, suggesting that the protein interior was more stable than the surface. Thus, we narrowed the selection of mutation sites to loop regions located on the surface of Anman.
There are many β-turns on loop regions, thus, we can design mutations based on the amino acid preferences on different residues (position i, i + 1, i + 2, i + 3) of β-turns (Supplementary Figure 3). In addition, proline has a higher frequency at position i + 1, it is the most rigid of all amino acids, and Xaa→Pro is beneficial for the improvement of thermostability Land et al., 2019). Thus, the mutation range was further narrowed to β-turns of Anman (types as furth, and we prefer mutation to proline. The Pymol software was used to calculate the average distance of C a for each β-turn from catalytic residues Glu168 and Glu276, and their correlation with the B-factor is shown in Figure 1. Except for β-turn 6, all other β-turns had a higher correlation in B-factor and distance from catalytic residues, and β-turns with a higher B-factor tended to be further away from catalytic residues. To obtain mutants that do not affect the catalytic efficiency, several regions (β-turn 3, 7, 9, 11) that were closer to the catalytic residues were excluded, and the remaining β-turns were used for the next round of screening.

Calculation of G and mutation
The remaining β-turns were subjected to virtual saturation mutation using the Discovery Studio 2019 software (Fu et al., 2017), and their G was calculated. The heatmap in Figure 2 shows the G for 380 mutants, and negative values represent mutations that are predicted to be stable. Since the calculation error for G was about 0.5 kJ/mol, a total of 315 mutants were predicted to be neutral or deleterious (their G > -0.5 kJ/mol). Some residues with high G values (including I12, Y154, R197, V300, and L312) were excluded. Finally, six single-point mutations (E15C, S65P, G83T, A84P, A195P, and T298P) were selected using a β-turn preference for flexible modifications. D114A and E222T (predicted to be negative mutations) were selected for use as the controls due to their higher 4 https://cpclab.uni-duesseldorf.de/cna/ G values (0.52 and 2.22 kJ/mol, respectively) and closer distances from catalytic residues (11.9 and 10.5 Å, respectively).

Construction of mutants and heterologous expression
Using the wild-type Anman gene as the template, all mutants were generated using polymerase chain reaction point mutation and verified via DNA sequencing. The primers used for all mutants are shown in Supplementary Table 1. Anman and positive mutants were overexpressed using P. pastoris X33 and then purified and verified by SDS-PAGE. The crude enzyme in lane 2 had an inconspicuous protein band at 80 kDa, while all of the remaining lanes showed a distinct single band at 47 kDa, which was inconsistent with the calculated molecular weights (37.1 kDa) due to the processing of glycosylation modification of Anman before its extracellular secretion by yeast cells (Supplementary Figure 4). Figure 3A shows the specific activity (black bars) and thermostability (red bars) of Anman and all mutants. Most of the mutants (excluding G83T and D114A) showed no significant decrease in a specific activity, which was due to the long distance between the mutation site and catalytic residue. As expected, the thermostability of mostly single-point mutants was improved, and A84P is the best-performing single-point mutant (from 46.2 to 68.7%). However, the relative activity of G83T was only 0.28-fold the level of Anman, which was an unexpected result. Finally, five single-site mutants with improved thermostability and no significant reduction in specific activity were combined to obtain the combinatorial mutant Mut5 (E15C/S65P/A84P/A195P/T298P). Mut5 had a 0.2-fold higher specific activity and a 0.7-fold higher thermostability than Anman (Table 1). It retained 82.4% of its activity after 10 min of incubation at 70 • C, which meant that the obtained thermophilic β-mannanase mutant had a significantly improved thermostability without the loss of enzyme activity. Heatmap of G for mutants. Green means G is negative and is predicted to be a stable mutation, while red means an unstable mutation, and the highlighted black squares indicate the mutation sites generated by this experiment.

Biochemical properties of Anman and mutants
To further investigate the performance of the mutant and Anman at high temperatures, their residual activity was measured after incubation at 70 • C for 0-60 min ( Figure 3B). The single-point mutant showed a slight increase over Anman, and E15C became the single-point mutant with the highest residual activity after 60 min of incubation, still retaining 22.3% of its specific activity, while Anman was almost completely inactivated. The combined mutant Mut5 had a remaining activity of 55.73%, which was significantly higher than that of Anman throughout the incubation process, indicating that Mut5 was more suitable for long-term function and industrial production at high temperatures than Anman.
Anman and all mutants showed the highest activity at 72.5 • C, with no significant change in the specific activity below 70 • C. But, the activity of Anman and mutants decreased dramatically when the temperature was increased to 75 • C ( Figure 3C). The mutation effect did not increase the optimal temperature of β-mannanase, which implies that it is more challenging to engineer thermophilic enzymes than mesophilic enzymes. The optimum pH of both Anman and the mutant was 3.0, which is close to the acidic environment of the animal stomach, suggesting that Anman and Mut5 can be used as feed additives in the breeding industry ( Figure 3D).

Multiple sequence alignment
To further explore the reasons for the decrease in enzyme activity or thermostability of G83T, D114A, and E222T, the protein sequence of Anman was aligned with 13 GH5 family thermostable β-mannanases, and their derived strains are shown in Supplementary  Table 2. The sequence conservation of E15, S65, A195, and T298 was less than 50%, G83 and E222 reached 92.3%, while D114 was up to 100% (Supplementary Figure 5). The comparison results indicated that G83, D114, and E222 were highly conserved residues, and their mutations occurred against the natural evolution, which was why they had negative mutation effects. The above results also demonstrated that relying on rational design alone could not completely exclude negative mutants, and introducing other methods in the screening process will be more beneficial in finding mutation sites. For example, Liu generated 44 mutants by modifying the surface-charged amino acids of ManAK, but only three mutants performed well .

Kinetic parameters and thermodynamic stability
To investigate the effect of mutation on the kinetics of Anman, the kinetic parameters of Anman and mutants were measured using 1-10 mg/ml LBG as the substrate (Table 1). Mostly, mutants had a similar K m as Anman, but E15C reached 6.7 mg/mL. Mut5 showed a slight increase in K cat , indicating that it had a higher affinity for LBG than Anman. There was no significant difference in catalytic efficiency K cat /K m for single-point mutants, but Mut5 was 0.25fold higher than Anman, reaching 239 mL/mg s (Supplementary Figure 6). The above results revealed that the mutation did not impair the catalytic ability of Anman for LBG, which further verified that Biochemical properties of Anman and Mut5. (A) Specific activity (black column) and relative activity (red column) after 10 min incubation at 70 • C of Anman and mutants, and relative activity was obtained by the ratio of residual activity after 10 min of incubation at 70 • C to the initial activity without incubation. (B) Relative activity of Anman and mutants after incubation at 70 • C for 0-60 min. (C) The optimal temperature of Anman and mutants. (D) Optimal pH for Anman and mutants. Relative activity † (%) mutations farther away from the catalytic residues had a lower effect on enzyme activity and that it could be one of the conditions used for screening thermostable mutants. The T m for Anman was 75.7 • C, which was close to its optimal temperature, and this was likely due to the protective effect of LBG on the catalytic pocket of β-mannanase. Except for A195P, the T m for other mutants was higher than that for Anman, but the T m for all single-point mutants was less than 1 • C, and the T m value for Mut5 was only improved by 2 • C compared to that for Anman (Supplementary Figure 7). Anman only had a 1.2% residual activity at 75 • C, while Mut5 retained 36.0% activity (Supplementary Figure 8), and both Anman and Mut5 were completely inactivated when the temperature reached 80 • C. This property made Mut5 more suitable for high-temperature conditions in industrial applications.
The t 1/2 was obtained by incubating at 70 • C for different amounts of time. The t 1/2 for A84P was the highest of all singlepoint mutants and 1.2-fold higher than that for Anman. Although the t 1/2 for S65P and A84P was significantly higher compared to that for Anman, Mut5 was 7.8-fold higher than that for Anman, reaching 73.4 min, and this phenomenon was caused not only by epistatic effects but also by the modification of most of the flexible sites of Mut5. When there are multiple flexible sites on the surface, a single-point mutation often cannot prevent the process of thermal inactivation; therefore, if multiple or all flexible sites on the surface can be rigidified, enzyme thermostability will be dramatically enhanced, which is why Mut5 performed much better than the single-point mutants.

CNA and MD simulation
Constraint network analysis can be used to visually compare the overall rigidity of Anman and Mut5. The P ∞ of Anman and Mut5 decreased sharply when E cut−hb was -1.84 and -2.76 kcal/mol, The constraint network analysis and molecular dynamics simulation of Anman (black) and Mut5 (red). (A) P ∞ changes during thermal simulation, and Anman and Mut5 change from rigid to flexible at E cut−hb of -1.84 and -2.76 kcal·mol -1 , respectively. (B) The number of hydrogen bonds of Anman and Mut5 during MD simulation. (C) RMSF of Anman and Mut5 at 370 K. Higher RMSF indicates more unstable residues, and the blue sticks represent the corresponding regions. Dynamics cross-correlation map of C α atoms in Anman. The white region represents residue motion without correlation; the purple region (cyan region) represents residue motion with correlation and in the same direction (opposite direction), and the darker color represents a stronger correlation.
Frontiers in Microbiology 07 frontiersin.org respectively ( Figure 4A). The simulated temperatures were 336.8 and 355.2 K, according to Equation (3), which were consistent with the previous experimental results for T m values. This also meant that Anman was inactivated earlier than Mut5. To further investigate the local and global stability of mannanase, a 50-ns kinetic simulation was carried out. It showed that the RMSD for Anman was usually higher than that for Mut5 ( Supplementary Figure 9), and Mut5 had more hydrogen bonds for the first 45 ns ( Figure 4B). These results indicated that the flexible engineered Mut5 was more stable than Anman in global conformation and more resistant to unfolding at high temperatures. MD simulation was also used to compare the local conformational stability of Anman and Mut5. The RMSF for all five mutation sites was reduced to various degrees, indicating that the mutation effect increased the rigidity of these sites ( Figure 4C). The RMSF values in the region near the P65, P84, and P195 mutation sites were substantially reduced because the pyrrolidine proline ring did not only make the C a and N atoms of the backbone more rigid in order to not rotate easily but also reduce the conformational flexibility of the nearby residues . In addition, regions 24-35, 54-58, 106-117, and 307-321 also significantly reduced RMSF in the absence of mutations, probably because the mutation effect caused rearrangement of the internal force network through long-distance perturbation pathways (El Harrar et al., 2021). To investigate this phenomenon, the dynamic cross-correlation matrices (DCCMs) were calculated for Anman ( Figure 5). There was a prominent cyan spot in the regions 54-58 and 106-117, which means that the atomic motions in these two regions were correlated and moved in the same direction, so their RMSF rise and fall trends should also be the same. In addition, region 54-58 was also spatially distributed adjacent to region 106-117, and there are many hydrogen bonds between the two regions; therefore, it is reasonable that the mutation causes a decrease in flexibility in both regions. The same process occurs in regions 24-35 and 307-321. The motion correlation of residues can be further used for the trade-off between thermostability and activity (Yu and Dalby, 2018).

Molecular force network analysis
The present study analyzed the chemical bonds within 5 Å of each mutation site ( Figure 6). Several hydrophobic residues were located around E15, and the hydrophobicity of cysteine was much higher than that of glutamate. Thus, the E15C mutant site formed hydrophobic interactions with the surrounding residues to stabilize the region. A similar phenomenon was observed in lipolytic enzymes (Gallardo et al., 2010). The N atom of C15 also formed an additional hydrogen bond with the O atom of I12 ( Figure 6A), which increased the thermostability of the E15C mutant. Even though A195P did not form more hydrogen bonds than Anman, D194 formed a salt bridge with R197 in A195P due to the structural perturbation caused by the mutation. To verify this phenomenon, the number of frames forming a salt bridge was divided by the total number of frames in the MD simulation. As a result, the probability of forming a D194-R197 salt bridge increased from 51.1 to 91.1%, and the distance between the OD2 atom of D194 and the NH2 atom of R197 decreased from 6.3 to 3.5 Å. Therefore, the improved thermostability of A195P was due to the formation of a salt bridge ( Figure 6B). For the T298P mutant, no hydrogen bonds formed since the angle between the O atom of T298 and the N and H atoms of V300 was lower than 120 • . In comparison, the mutated proline changed the local conformation and formed hydrogen bonds between P298 and V300, thereby reducing the flexibility of the local region ( Figure 6C).
There were no hydrogen bonds near the S65 region of Anman, and there were four hydrogen bonds in S65P. The O atom of S64 formed two hydrogen bonds with the N and OG1 atoms of T67 ( Figure 6D), which made the unfolding of the S65P region require more heat. Therefore, S65P did not easily deactivate at high temperatures. The same happens in the A84P region ( Figures 6E, F), where the NE2 atom of Q88 approached the T82 residue to form an additional hydrogen bond, and the OD1 atom of D85 and the O atom of P84 formed two hydrogen bonds. A situation where one atom becomes an acceptor for two hydrogen bonds can strengthen the binding between residues. In addition, the S65P and A84P mutation sites were relatively close to each other in space, and the mutation effects could be transmitted between them (Dinmukhamed et al., 2020). Therefore, the force network between the S65P and A84P regions could be significantly changed or even redistributed, and the mutational effect resulting from structural rearrangement may also have a positive epistatic effect, which may account for the significantly reduced RMSF of residues near the two regions.

Conclusion
In conclusion, β-mannanase derived from A. niger CBS513.88 was characterized and its biochemical properties were measured. The combined mutant Mut5 had an apparent increase in t 1/2 , which was 7.8-fold higher than that of Anman. Meanwhile, the T m of Mut5 was also increased by 2while, the ed. The combined mutant Mut5 had to determine that Mut5 had more global and local stability. Additional chemical bonds were generated near the corresponding mutation sites, which is also the molecular mechanism behind the better Mut5 thermostability compared to that of Anman. The present study yielded a β-mannanase mutant that was more suitable for industrial applications. The RFS strategy also provided an effective way for screening enzyme hotspots.

Data availability statement
The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.